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"Characterizi ng surfaces in medical imaging" 
Description 

Field of the Invention 

The invention relates to an image processing system having image data 

5 processing means for characterizing boundaries of regions in an image and for matching 
surfaces in medical imaging. The invention also relates to a medical examination 
apparatus having such an image processing system. The invention further relates to an 
image processing method having steps for operating the system and to a computer 
program product having instructions for carrying out the method steps. The invention 

10 finds its application in the field of medical imaging and, more especially, in the field of 
x-ray medical imaging. 

When analyzing 3D medical images, it is often useful to match a discrete surface 
model to the boundary of an anatomical object. This can help, for example, in extracting 
the object from other unrelated objects for improving the visualization of said object. 

15 Among many other applications, such discrete surface model allows to make clinical 
measurements such as linear dimensions, boundary surface area, cross-sectional area, 
enclosed volume of an organ. It is also often desirable to compare or match active 
surface models with respect to one another. 

Background of the Invention 

20 A publication entitled "Edge Dipole and Edge Field for Boundary detection", 

pp. 179- 189, in "Part of the SPIE Conference on Hybrid Image and Signal Processing 
VI, ORLANDO, Florida, April 1998, SPIE Vol. 3389", by Toshiro KUBOTA et aL, 
proposes a framework for treating edges as directional dipoles that induce a field around 
themselves. This paper treats edges not only as vectors but also as dipoles, which induce 

25 a vector field around them. Thus, said dipoles are called edge dipoles and the field is 
called edge field. An analogy is made between this concept and the interaction of 
magnetic dipoles with the magnetic field. The field induced by an edge dipole exhibits 
smooth continuous flows from the positive pole to the negative pole. The field in turn 
can be used to rotate other dipoles to align with the field. The approach of this paper 

30 provides an interaction mechanism between the dipole and the field, which has an 
influence on the edges. Also, the edge dipole field has circular flows diverging at the 
positive pole and converging at he negative pole, which provides a convenient 
mechanism for lateral inhibition and discourages thick boundaries from forming. The 
dipoles interact with the field and align themselves into a smooth contour configuration. 
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This concept is used in edge linking in order to complete missing image information. 
So, according to the disclosed method, it is needed to compute the influence of the 
dipoles within the region surrounding them in the images themselves. 

Summary of the Invention 
5 The present invention has for an object to provide an image processing system 

having processing means for characterizing the boundaries of regions using edge 
dipoles without computation related to the image region itself. The present invention 
has further for an object to provide such an image processing system for analyzing an 
anatomical surface of interest for comparing, matching or adapting a geometric model 
10 with anatomical object boundaries observable in a medical image, and/or a geometric 
model with another geometric model, and/or real image data with other image data of 
anatomical object boundaries observable in a medical image. Using the proposed 
system, substantial portions of the objects can be matched together in a patch wise 
approach. The output of the present system is valid whether the objects are virtual, such 
15 as active models, or are real such as represented by real image data. A basic principle of 
the invention is to represent a boundary of an object by dipoles alike a set of small 
elemental magnets. Each dipole characterizes an element of the boundary. Each dipole 
has a strength proportional to the size of the boundary element, is normal to said 
boundary element and oriented toward a direction regarded as outward direction with 
20 respect to the object limited by the boundary. This system means are recited in Claim 1 . 

The image processing system can be implemented as a specially programmed 
general-purpose computer. The image processing system can be a workstation. 
According to the invention, an image processing method having steps for operating such 
this system is also proposed. The present invention yet further provides a computer 
25 program product having a set of instructions, when in use on a general-purpose 
computer, to cause the computer to perform the steps of this method. The present 
invention still further provides a medical examination apparatus incorporating the image 
processing system putting into practice this method in order to process medical image 
data obtained by the imaging apparatus, and means for visualizing the image data 
30 produced by said method. The visualization means typically consists of a monitor 

connected to the image processing system. Advantageously, the workstation and image 
processing system of the present invention are interactive, allowing the user to influence 
clinical data that are evaluated and/or the manner in which evaluated data is to be 
visualized. 
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Brief Description of the Drawings 

The invention and additional features, which may be optionally used to 
implement the invention to advantage, are described hereafter with reference to the 
schematic figures, where: 
5 Fig. 1 A illustrates dipoles associated to a boundary of an object and FIG. IB 

illustrates the definition of the polynomial function parameters; FIG.1C illustrates a 
Dirac function; 

FIG.2 is a block diagram illustrating elements of a medical examination 
apparatus, incorporating a medical viewing system; 
10 Fig.3 A is a flow diagram showing the main steps of an image data processing 

method for characterizing a boundary and FIG.3B is a flow diagram showing the main 
steps of an image data processing method for matching two boundaries; 

FIG.4A illustrates a starting mesh model and FIG.4B illustrates a step of fitting 
the model to a heart cavity; 
15 FIG-5 illustrates the mesh model that matches the heart cavity. 

Description of Embodiments 
An image processing method for characterizing boundaries of regions in an 
image is 

first described. The invention further relates to the application of this method to process 
20 medical images for improving the visualization of an anatomical surface of interest. The 
present image processing method may be used for analyzing an anatomical surface of 
interest for comparing, matching or adapting a geometric model with anatomical object 
boundaries observable in a medical image, and/or a geometric model with another 
geometric model, and/or real image data with other image data of anatomical object 
25 boundaries observable in a medical image. Matching together substantial portions of the 
objects can be performed in a patch wise approach. The present invention preferably 
makes use of mathematical tools based on a polynomial function such as the Hermite 
transform, also known as Hermite wavelets. This mathematical tool is succinctly 
described hereafter to facilitate computation steps that are advantageous to carry out the 
30 method of the invention. 

In a 2D or 3D image, the image intensity f is defined in function of the location 
of each image pixel or image voxel. So, in a 2D image, the scalar value image intensity 
f is given as function of the real positional coordinates x,y . A positive valued function 
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w(x, y) is introduced to define a fuzzy observation window W in the 2-D image. 
Geometric moments are defined as a set of quantities: M nvn 

= JJw(x, y) f(x, y) x m y n dx dy (1) 

The whole family of moments provides an alternative representation of the image 

5 f(x, y) that fully and uniquely characterizes the intensity fi[x, y) within the observation 
window w(x,y) . In practice, lower orders of the data (m,n) are the most robust 
whereas higher orders are more delicate to compute numerically and are more noise- 
sensitive. These moments are powerful in coding morphological information, but they 
have two shortcomings: 

10 a) When comparing two images together, within a similar observation window, 

it is not obvious how similarity/ dissimilarity measures of their moment representation 
relate to similarity/ dissimilarity of the original images. 

b) It is tedious to work out a comprehensive set of related measures that are 
invariant under a transformation group. 

15 In order to solve these shortcomings a) and b), instead of using monomials x'V* 

in (1) one can use polynomials P^ n (x, y) , which are mutually orthogonal within the 
observation window w(x, y) in the sense that: 

= JHx,y)P n , n (x,j;)P lrf . rf (x,y)dx.dy (2) 
where each 5^ and are equal to 1 if the two integers indices (m,iri) or (n,n') 

20 coincide and to 0 otherwise, and where > 0 is referred to as the norm of the 
polynomial P^fo y) . In fact it is possible to build a complete set of orthogonal 
polynomials provided the window function w(x,y) is never negative. Such a 
polynomial basis can provide least mean square approximation of an arbitrary function 
f(x,y) in the form: 

25 f(x,y) = S f n, n ^(x.yVC^ (3) 

m,n 

where the 2-D Hermite coefficients f^ n are given by: 

f n,n - JJ^x f y).P^(x,y).l[x t y).dxdy (4a) 
Here, the coefficients f^ n replace the geometric moments M nvn defined in the previous 
paragraph. Those new coefficients are referred to as orthogonal moments. An important 



WO 2004/023396 PCT/IB2003/003775 

property answers the requirement a) stated at the beginning of the current section. This 
theorem, called ParsevaFs theorem, says that the scalar products are conserved. So, 
considering any two real-valued images f(x,y) and g(x,y) with respective orthogonal 
moments f^ n and g^ , it is seen that taking the scalar products in transformed space, 
5 where the orders m and n play the role of coordinates, is the same as taking the product 
in x,y space: 

IXn g,,^,, = JJw(x,^)f(x,y)g(x,y)dxdy (5) 

m,n 

This relationship can be deduced from the orthogonality of the polynomial functions 
P^ n . It is noteworthy that the normalization factors 1/C^ play in the m,n space, the 
10 role that the window function w(x,y) plays in the x,y space. 

The second shortcoming b) of the geometric moments is that the moments do 
not transform easily with changes of referential. Therefore, it is required to simplify 
transformations of the orthogonal moments with change of referential. 

An advantageous way of doing this is to choose Hermite polynomials. This 
15 corresponds to taking for the window an isotropic Gaussian function centered about any 
reference point (£,77) , so that: 

This resulting transform is called the Hermite transform. The corresponding 2-D 
polynomial basis can be expressed in term of products: 
20 P ra , n (x,y) = H m (x).H n (y) (7) 

where H m and H n are Hermite polynomials of orders m and n respectively. The 
corresponding normalization factors are: = 2 m+n m!n! 
(8) 

Since an object of the invention is to compare images together, it is needed to 
25 easily deal with geometric transformations. 

Dealing with translations is first described: The Hermite coefficients f m>n of 
function fQ can be expressed in terms of partial derivatives of a smoothed function F() 
obtained by convolving fQ with the Gaussian kernel w(). Those skilled in the art can 
easily demonstrate that: 
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m+n ps 



dx m dy' 



(9) 



This form allows to obtain the effect of translation by a simple application of 
Taylor theorem where the coefficients are directly related to the Hermite transform 
coefficients f m>n . The fact that a Gaussian smoothing was obtained to evaluate F implies 

5 that it is enough to use very few terms for the Taylor expansions provided one stays 
within a disk of radius 1 or 2a from the original window centre. 

Dealing with rotations is then described: To deal with 2D rotations, the most 
suitable related transform is that of Gauss-Laguerre transform described in the 
following publication by G. Jacovitti and A. Neri; "Multiresolution Circular Harmonic 

10 Decomposition", IEEE Transactions on Signal Processing, vol. 48, pp. 3242-3247 
(2000). This transform may be related to that of Hermite since the corresponding 
polynomials are related to solutions of the Quantum harmonic oscillator as described in 
the following publication by J. J. Koenderink and A. van Doom; "Generic 
Neighbourhood Operators'*, IEEE Transaction on Pattern Analysis and Machine 

15 Intelligence, vol. PAMI-14, pp. 597-605 (1992). Hermite polynomials result from 
seeking solutions using Cartesian coordinates whereas the Laguerre polynomials 
appearing in the Gauss-Laguerre transform result from the use of polar coordinates. The 
essential thing one needs to know is that, for any given order (psm+n) the two sets of 
solutions can be transformed into each other. The Cartesian set of polynomials can be 

20 written as products Pnwi(x,y)=H m (x).H n (y) whereas the polar form (Gauss-Laguerre) are 
in the form 

^ V '[sin(*0)J 



where Q p ,k(r) is a function of the radial coordinate r = <yjx 2 + y 2 and 9 is the angular 
polar coordinate. The order k of the circular harmonic factor spans all integers in the 

25 interval [0, p] with same parity as p (for k=0 only one circular harmonic exists, namely 
the constant 1). It is easily seen that, for both representations, there are exactly (p+1) 
independent functions. In order to find the (p+1) 2 matrix that transforms the Hermite 
into Gauss-Laguerre functions, it is enough to compare angular behaviours for very 
large r values. Elementary trigonometry can be used to express P m , n (x,y) °c 

30 cos m (9).sin n (0) in terms of the circular harmonic functions referred to earlier. It is also 
easily established that the same matrix can be used to transform coefficients from one 
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basis to the other. One can also deduce the new normalization factors Cm^ from this 
matrix because both set of functions form orthogonal bases. In order to revert from 
Gauss-Laguerre back to Hermite, the inverse matrix is used. 

Rotation with an angle A0 is easily implemented because the coefficients C p ,k 
and Sp, k related to the Gauss-Laguerre functions Q Pi k(r).cos(k6) and Q p ,k(r).sin(k0) 
behave exactly like x and y components of a 2-D vector rotating with an angle kA9 (for 
non-zero angular harmonic order k). Thus, a rotation would transform C Pt k and S p ,k to 
CVk and S' p ,k such that 

cos(£A0) sin(*A0)YC v 
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Which is totally equivalent to the complex variables equation: 

(C;, +/5;,J=exp(i^A^)(C^ +iS pJt ) (12) 
where i stands for the imaginary unit and A9 stands for the rotation angle of the 
referential. 

Dealing with scale changes is further described: Changing scale implies 
changing the Gaussian window function. However, one can alternatively imply the scale 
changes in the approximating polynomials P m , n 0- Although this approach is relatively 
simple, it implies using a least-mean square approximation obtained within a given 
window for extrapolation into a larger window or for interpolation within a smaller 
window. Scale changes are most easily operated on monomials. The approach proposed 
is just to express the polynomials as a sum of monomials, operating scale change and 
then reverting back to the orthogonal representation. For Hermite polynomials, this can 
be expressed as matrix operations in the form H' = AhA~ x H , where: 
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(13) 



and where A is the lower triangular matrix allowing to express the column vector of 
Hermite polynomials in terms of the column vector of the powers of x. This allows to 
express H n (A.x) as a linear combination of original Hermite polynomials H p (x) of equal 
or smaller order (and of same parity). The matrix A A A' 1 can be computed once for all, 
it is a sparse lower triangular matrix in which all non-zero coefficients are polynomials 
in X of orders never exceeding the order of the line. Since 2-D Hermite polynomials are 
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just products of ID polynomials, this scaling operation is dimensionally separable. 
Exactly the same transformation applies to Hermite coefficients. Because of the 
sparseness of the matrices involved and of dimensional separability, scale changes are 
not computationally costly. 

5 The above-described operations can be applied in images of more than two- 

dimensions: Most of the above-described operations easily generalize to any number of 
image dimensions. The only basic change is how one deals with rotations. Spherical 
harmonics have to replace circular harmonics. For 3-D, the corresponding spherical 
harmonic formulae can be found in classical mathematics handbooks and the conversion 

10 matrices from Cartesian (Hermite) to polar form (Laguerre) can recursively be 

evaluated as is done in 2-D. However, a rotation of the referential implies the use of 
difficult spherical harmonics addition formulae. Another formulation has been disclosed 
by the Clifford Algebra group of Ghent University by R. Delanghe, F. Sommen, V. 
Soucek, F. Brackx in "Clifford Algebra and Spinor- Valued Functions, A Function 

15 Theory of the Dirac Operator", Kluwer Academic Publishers (1992). For the 3-D case, 
the approach proposed by this publication just needs the introduction of unit called 
"quaternions" to implement 3-D rotations of the Gauss-Laguerre coefficients, thus 
providing a generalization of the complex rotation factor cxp(ikA0) as discussed in 2D 
case. Now, quaternions are entering the engineering field and are not considered to be 

20 exotic anymore in the image processing community, as described by B. Horn in "Closed 
Form Solutions of Absolute Orientation using Unit Quaternions" in Journal of the 
Optical Society, vol. A4, pp. 639-642 (April, 1987). The basic advantage in using the 
above set of proposals is to transform a set of image processing operations into 
hopefully simple algebraic forms. These mathematical tools may be used by those 

25 skilled in the art in order to perform the computation steps of the image processing 
method of the invention. 

A method that applies this technique of characterizing a boundary using the 
function of intensity with the Hermite coefficients is described hereafter. The present 
image processing method has steps for computing Hermite coefficients in order to 

30 characterize a curve line or a boundary of an object in an image window by the Hermite 
transform. This method permits of modelling and extracting interfaces or boundaries of 
objects in medical images. For performing this operation, it is possible to compute 
Hermite Transform Coefficients of an object boundary regarded as an oriented interface 
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by sampling the corresponding boundary. There is no need to get region data as 
described in the cited "KUBOTA" publication. 

FIG.1A represents schematically a region of interest ROI of an object in a two- 
dimensional (2-D) image. The region of interest ROI shows a 2-D boundary. A window 

5 W defines a portion SOI of the boundary that is considered. The 2-D boundary SOI is 
decomposed into small adjacent elementary linear elements, called segments, such as 
SI, S2, S3, S4, S5, respectively related to a predetermined number of boundary pixels, 
called reference pixels. Every segment such as S1-S5 of the boundary SOI in a 2D 
image is represented by a respective dipole Dl, D2, D3, D4, D5. The arrows Dl to D5 

10 illustrate the boundary dipoles. In 3-D images, the region of interest ROI shows a 3-D 
surface of interest SOI, called 3-D boundary, which is decomposed into small portions 
in a window W. Every small portion is represented by a dipole of strength proportional 
to the area portion of the 3-D boundary and oriented along the outward normal 
direction, at a defined centre of the small portion. 

15 From such a representation, the Hermite transform can be computed within the 

window W as schematically represented in FIG.l A for a 2-D image. The intensity of the 
boundary SOI in the window W is defined by a function of intensity f(x,y) where x,y 
are the coordinates of the reference pixels at the center of the segments or small 
portions. Assuming that the window centre is at the origin and that distances are 

20 normalized so that <W2=7, the Hermite transform coefficients related to a dipole vector 
(I x ,I y ) at point (x,y) can be expressed as: 

fm,n = w(x,y).(I x .H m +i(x).H n (y) + I y .H m (x).H n +l(y)) 

(4b) 

The Hermite transform of the interface dipole layer can be approximated as the sum of 
25 such contributions over all discrete interface elements in the window. Generalization to 
more than 2-D is again straightforward. 

A first data, denoted by n, called spatial resolution, is related to the number of 
said segments in which the boundary is decomposed in said window W. A second data, 
denoted by m, called intensity resolution, is related to the number of possible gray 
30 levels for the pixels of the boundary in said window W. The definition of the data n, m 
is illustrated by FIG. IB. Near point O, the spatial resolution n is low, meaning the 
number of segments is small while the segments are long; and the intensity resolution is 
low, meaning the number of gray levels is small. Instead near point p, the spatial 
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resolution n is high, meaning the number of segments is large while the segments are 
small; and the intensity resolution is high, meaning the number of gray levels is large. 
The data n, m may be chosen a priori by the user. So, using the values of n, m near 
point O, the boundary SOI in the window W is characterized by a few number of long 

5 segments and a few number of gray levels. Instead, using the values of n, m near point 
p, the boundary SOI in the window W is characterized by a large number of small 
segments and a large number of gray levels. 

In the present method, instead of regarding a boundary pixel directly, said 
boundary pixel is first coded using the above described technique. Each segment is 

10 characterized by a dipole, such as Dl, D2, D3, D4, D5, located at the centre x,y of the 
segment, oriented along the outward normal direction of the object and having a 
strength that is defined using the spatial and intensity data related to the segment. The 
coding data include the first information related to the spatial resolution n and the 
second information related to the intensity resolution in the window W of the region of 

15 interest ROL The spatial resolution n is considered as a first dimension and the intensity 
resolution m is considered as a second dimension for determining coding data related to 
a pixel at the location x,y . Alternately as shown in FIG. 1C, the segment may be coded 
by a Dirac intensity function having a highly positive value, outside the object, with 
respect to an axis directed along the dipole, and a highly negative value inside the 

20 object. 

The major difference between the present invention and the technique disclosed 
in the publication "KUBOTA et aF cited above, is that instead of broadcasting the 
effect of each dipole in a portion of the image as was done in the known technique, the 
effect of each such dipole is integrated into transformed coefficients. The effects of 

25 those dipoles, such as D1-D5, within the window W, are cumulated as coefficients of an 
image transform procedure. The polynomial Hermite transform is advantageously used. 
As above-described, said Hermite transform is capable of coding the boundary inside 
the window using a finite number of coefficients. The contribution of each dipole 
element can be computed analytically knowing its position within the window, its 

30 strength and its orientation. The coding data form a series of two-dimensional (2-D) 
coefficients for defining said polynomial function. The number of 2-D coefficients is 
small when the coding data are chosen near point O of FIG. IB, while the number of 2- 
D coefficients is large when the coding data are chosen near point p of FIG.1B. The 
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number of 2-D coefficients may be chosen within a range in which it may vary slowly, 
since the coding data n, m may be each chosen in a range in which they may vary 
slowly. According to the method of the invention, the boundary is characterized by said 
polynomial function with said 2-D coefficients in the window W. If the polynomial 

5 function is constructed with a great number of such 2-D coefficients, the representation 
of the boundary is smooth. If the polynomial function is constructed with a small 
number of coefficients, the representation of the boundary is coarse. Hence said 
polynomial function with said 2-D coefficients permits of classifying the representation 
of the boundary in a hierarchical manner, from coarse to fine, according to the number 

10 of coefficients used to construct the polynomial function. The medical viewing system 
and an image processing method of the present invention permit to minimize the load of 
computation produced by the KUBOTA et al approach. 

FIG.3A is a flow diagram of main steps of the present method, comprising: 

In step SO, acquisition of image data: The image data input to the method can be 

15 for example 3-D computed tomography image data obtained for a subject heart. The 
medical image data consists of a large number of data relating to points, each 
corresponding to a respective position within the patient's body. The image data can be 
submitted to preprocessing in order to eliminate noise. The method further comprises 
steps of: 

20 In step S 1 , computation of image data of an object surface denoted by B 1 or B2. 

In step SI, the outer surface of the heart muscle, denoted by Bl, is identified from 
within the image data via a segmentation process as illustrated by the segmented 2-D 
curved line or 3-D surface SOI in Fig.lA. In another technique, the 3-D surface, 
denoted by B2, may be obtained as an active model providing a best fit to the heart 

25 muscle, or other anatomical object under consideration. Techniques for producing an 
active model of an anatomical object are well-known, for example by the description in 
the publication entitled "General Object Reconstruction Based on Simplex meshes" by 
Herve Delingette, in the International Journal of Computer Vision, 32, 1 1 1-142, 1999. 
Such a model is illustrated by FIG.4A and FIG.4C. 

30 In step S2, decomposition of the boundary into 2-D elements of a 2-D curve line 

or small portions of a 3D boundary, such as Bl or B2, as above described in reference 
to FIG.1A. 
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In step S3, characterization of the 2-D or 3-D boundary elements for example by 
dipoles or Dirac functions as illustrated by FIGs.l A-1C; and computation of 
corresponding Hermite coefficients according to formulae (4a) or (4b). 

In step S4, cumulating the effects of the dipoles Dl to D5 within the window W 

5 as coefficients of an image transform procedure. For example, one can use a Hermite 
transform that is capable of coding an image inside the window W using a finite number 
of coefficients. In this step, the contribution of each dipole element Dl to D5 in the 
window W is computed analytically knowing its position within the window, its 
strength and its orientation according to formulae (3), (4), (6) and (7). 

10 In step S5, computation of the Hermite coefficients for translation, rotation and 

scale cange according to formulae (9), (1 1), (12), (13). 

The Hermite transform of a boundary can now be used for different purposes. 
The present invention proposes a method having the steps described above and having 
further steps for matching two surfaces. For example this method has further steps to 

15 compare and/or to match a discrete model interface to another discrete interface; or to 
compare and/or to match a discrete interface to a real image. It is noted that the above 
proposals can substantially improve the performance and computational complexity of 
matching procedures. None of the above procedures requires to know a mapping from 
initial point position (prior to transformation) to a target position (after transformations) 

20 as is needed by the cited algorithms disclosed in the publication by B. Horn. So, the 
invention further relates to applying the first steps of characterizing the function of 
intensity of a boundary with the Hermite coefficients for comparing or matching two 
objects defined by curve lines or surfaces in corresponding windows. The operation of 
comparing or matching the boundaries of said two objects can be done only using the 

25 respective Hermite coefficients relating to their boundaries. The transformation needed 
to match a portion of interface with another one (registration or motion estimation) can 
be estimated for a pair of objects of same nature such as model with model or real data 
with real data; or of dissimilar nature such as model to data or data to model. The 
proposed method can be applied to medical image processing in order to enabling 

30 improved matching of surfaces such as an active model with anatomical object 

boundaries observable in a medical image, and/or an active model with another active 
model, and/or real image data with other image data of anatomical object boundaries 
observable in a medical image. 
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Hence, this method comprises steps for computing Hermite coefficients in order 
to characterize corresponding boundaries of two objects in an image window by the 
Hermite transform, including Hermite coefficients for translation, rotation, and scale 
changing and further steps for comparing or matching said two boundaries. Here, the 

5 template (or block-) matching, which is well known to those skilled in the art, can be 
transformed into matching of transformed coefficients. Rather than making exhaustive 
searches for translations, rotations, and scale changes in said template (or block-) 
matching operations, the mean square matching error can be written as an algebraic 
function of the transformation parameters. Optimal match between two boundaries can 

10 then be found by conventional gradient descent. Furthermore, it can be hoped that the 
simple algebraic form of the matching error can be exploited to efficiently locate all the 
local minima. A rather interesting feature of the Hermite transform is that the set of 
coefficients for which one of the indices is 0 corresponds to the Hermite Transform of a 
"Weighted Radon" projection along axes corresponding to the 0 index. This may be 

15 useful, for example, for matching a 3-D pattern with a 2-D projection pattern. 

Referring to Fig.3B that is a flow diagram, this method comprises: 
In steps Tl, T2: respectively acquiring image data of first and second images. In 
the example described hereafter, in the first image the object is an organ represented by 
real data, and in the second image, the object is a virtual model of the organ having a 

20 boundary to be compared with a boundary of said organ. The interface of the virtual 
model is a binary region. 

In steps T3, T4: respectively segmenting said images to provide image data of 
the corresponding boundaries respectively denoted by Bl, B2 in the first and second 
images. 

25 In step T5 : calculating Hermite coefficients of the Hermite transforms for 

boundaries Bl and B2, as described in reference with steps S2, S3, S4 illustrated by 
FIG.3A. 

In step T6: deducing, as in step S5 illustrated by FIG.3A, coefficients resulting 
from rotation, scale-change or translation of one object with respect to the other. The 
30 computations only require knowledge of the geometric transformation and of the 
Hermite coefficients before transformation. 

In step T7: estimating the transformation needed to match the portion of 
boundary or interface of the model with the portion of boundary or interface of the 
organ, only using their respective Hermite coefficients. 
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Another application is for performing registration or for motion estimation. 
Another application consists in smoothing an interface model through adapting it to a 
flat, or another appropriately shape target interface. FIGs.4A illustrates an active model 
to be matched to heart cavity as shown in FIG.4B. This model is deformed using 

5 internal and external forces in such a way as to fit the heart cavity. In order to determine 
the best match, the model surface boundary, as shown in FIG.5, is compared to a 
segmented surface boundary of the heart cavity using the above-described method. In 
this first example, the above-described steps are used for estimating the transformation 
needed to match a portion of interface of objects of dissimilar nature, such as a model 

10 compared to real data or such as real data compared to a model. The above-described 
steps can also be used for estimating the transformation needed to match a portion of 
interface of a pair of objects of same nature, such as a model compared to a model or 
such as real data compared with real data. Then, each reference point of the reference 
surface is processed in turn and the normal to the reference surface is calculated at each 

15 reference point. All operations involved can be done in computationally efficient 

manner. The basic technique proposed here is generic; more than the cited operations 
may be carried out. 

In general, the described method above, which is based on the technique of 
characterizing the function of intensity of a boundary with the Hermite coefficients of 

20 the Hermite transform of the boundary, can be used for purposes such as: 

Adapting one contour to an image, which would consist in seeking a maximum 
of the correlation. This operation adapts the position of the dipoles in such a way as to 
maximize correlation of an interface with another interface called target interface. This 
can be used for example to model large portions of the boundary of an anatomical 

25 object using an active contour or active surface as model. This can be done by changing 
the dipole positions within the window (free-form adaptation) or by seeking optimal 
rotation, scale change and/or translation (global adaptation within the window). In fact, 
the operation of comparing or correlating the boundaries of said two objects can be 
done only using the respective Hermite coefficients relating to their boundaries. 

30 Smoothing an interface model: this is possible through adapting it to a flat or 

appropriately shaped target interface. Smoothing such a boundary is performed by 
minimising the amplitude of higher order terms of the Hermite transform. 

In another application, this method can be used to study noise or texture in 
images. This operation can be carried out using estimation of their autocorrelation. This 
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consists in studying how the image correlates to a translated version of itself. An 
extension of this idea is to characterize the way the image correlates to itself when it 
undergoes, translation, scale changes and rotations. One should therefore estimate: 

R(T) = [f(T~ x x)f(Xx^j 9 where T and T 1 are a transformation and its inverse. 

5 The windowed correlation product can be expressed algebraically in terms of the 

parameters defining the group of the geometric transformations T. If this group includes 
rotation, scaling and translation, the resulting autocorrelation function is rotation, scale 
and translation independent. If this is expressed in a functional form, coefficients 
characterizing this functional form can be extracted. For rectangular windows and when 

10 limiting T to mere translations, the functional form of the autocorrelation is 

characterized by the power spectral distribution (square of moduli of the Fourier 
transform) through the Wiener Khintchine theorem as described by W.H. Press, S.A. 
Teukolsky, W.T. Vetterling and B.P. Flannery; "Numerical Recipes in C, The Art of 
Scientific Computing", Cambridge University Press, 1992 (2 nd Edition). Such 

15 characterisations of autocorrelation functions can be used to characterizing texture or 
the spectral behaviour of noise, and or providing an invariant set of characteristics of 
shapes for pattern recognition. 

The above method for coding boundaries can be extended to other structures and 
to many applications. For example, oriented linear or tubular structures can be coded as 

20 a set of quadripoles. It is possible to put this transform in one of two forms: Hermite 

Cartesian form that is best suited to deal with translation or scale-change, and/or Gauss- 
Laguerre Polar form which is best suited to deal with rotations. Conversion from one 
form to the other is not computationally costly if one concentrates on low orders (i.e. 
when neglecting high frequency noise-prone terms in the expansions). Implementing 

25 these basic transformation operations can be done in a computationally efficient manner 
and can also be put in the form of compact algebraic expressions. Many practical 
applications can be envisaged as outlined in the previous section for putting the 3D 
rotation transformations into a mathematically correct framework and or characterizing 
the Generalized Autocorrelation function and relating it with invariant theories used in 

30 pattern recognition. An important advantage of Hermite transform is the ease and 
efficiency with which it can be put in forms that are particularly suited to deal with 
translation, rotation and scale change. Once the coefficients of the Hermite transform 
are known, one can deduce the coefficients resulting from rotation, scale-change or 
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translation of the object. The computations only require knowledge of the geometric 
transformation and of the Hermite coefficients before transformation. 

Fig.2 shows the basic components of an embodiment of an image processing 
system in accordance to the present invention, incorporated in a medical examination 

5 apparatus. As indicated schematically in Fig.2, the medical examination apparatus 
typically includes a bed 10 on which the patient lies or another element for localizing 
the patient relative to the imaging apparatus. The medical imaging apparatus may be a 
CT scanner 20. The image data produced by the CT scanner 20 is fed to data 
processing means 30, such as a general-purpose computer, that carries out the steps of 

10 the method. The data processing means 30 is typically associated with a visualization 
device, such as a monitor 40, and an input device 50, such as a keyboard, pointing 
device, etc. operative by the user so that he can interact with the system. The elements 
10-50 constitute a medical examination apparatus according to the invention. The 
elements 30-50 constitute a image processing system according to the invention. The 

15 data processing device 30 is programmed to implement a method of processing medical 
image data according to invention. In particular, the data processing device 30 has 
computing means and memory means to perform the steps of the method. A computer 
program product having pre-programmed instructions to carry out the method may also 
be implemented. 

20 The present invention is applicable regardless of the medical imaging 

technology that is used to generate the initial data. For example, when seeking to 
visualize the heart, magnetic resonance (MR) coronary angiography may be used to 
generate 3D medical image data in a non-invasive manner. See, for example, "Non- 
invasive Coronary Angiography by Contrast-Enhanced Electron Beam Computed 

25 Tomography" by Achenbach et al, in Clinical Cardiology, 21 , 323-330, 1998. The 

Achenbach et al article includes useful information regarding optional data processing 
steps that can be applied to the medical image data, for example, segmentation to enable 
a representation of certain anatomical features in isolation from others, etc. These steps 
can be applied in the method of the present invention. The present invention is 

30 applicable regardless of the way in which a surface of interest is modeled, whether via 
use of a reference simplex mesh, or in some other way. Various modifications can be 
made to the order in which processing steps are performed in the above-described 
specific embodiment. The above-described processing steps applied to medical image 
data can advantageously be combined with various other known 
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processing/visualization techniques. The drawings and their description hereinbefore 
illustrate rather than limit the invention. It will be evident that there are numerous 
alternatives that fall within the scope of the appended claims. Moreover, although the 
present invention has been described in terms of generating image data for display, the 
5 present invention is intended to cover substantially any form of visualization of the 
image data including, but not limited to, display on a display device, and printing. Any 
reference sign in a claim should not be construed as limiting the claim. 



